#MT Analysis of Survey Disposition
#PORES AAPOR
#December 2020

rm(list=ls())
#Set to the folder you've downloaded the material to
setwd()
load(file="DispositionDataForAnalysis.Rdata")
library(lfe)
library(stargazer)
library(survey)



#########
#Figure 1
#########

m1 <- felm(coop ~  log(cases) + dem.perc16 + density +
             female  +Age3039 + Age4049 + Age5059 + Age6074 + Age75p +
             other.race + black + hispanic + likely.party | state | 0 | state, data=dat)

coop.coefs <- coef(m1)
coop.ses <- m1$se
ypos.line <- seq(2,28,2)
ypos.coop <- seq(1,30,2)

labels <- c("ln(COVID Cases)","Dem % Pres. Election 2016",
            "Suburban Zip","Urban Zip","Female",
            "30-39","40-49","50-59","60-74","75 Plus","Other Race",
            "Black","Hispanic","Independent/Other","Republican")


pdf(file="./Figures/CLINTON-LAPINSKI-TRUSSLER-Figure1.pdf")
par(mar=c(5.1,12,4.1,2.1))
plot(coop.coefs, ypos.coop, xlim=c(-.1,.1), ylim=c(0,31),pch=16, col="darkorange", type="n", axes=F,
     ylab="",xlab="Effect on P(Cooperate)")
axis(side=2, at=ypos.coop, labels=labels, las=2)
axis(side=1, )
abline(v=0, lty=3, lwd=2)
abline(h=ypos.line, lty=2)
points(coop.coefs, ypos.coop, pch=17)
segments(coop.coefs, ypos.coop, coop.coefs+1.96*coop.ses, ypos.coop)
segments(coop.coefs, ypos.coop, coop.coefs-1.96*coop.ses, ypos.coop)
dev.off()

#Substantive impact in PA

#How many Republicans in PA sample?
sum(dat$coop[dat$state=="PA" & dat$likely.party=="R"],na.rm=T)

#How many Republicans if participation rate rose to level of Democrats?
(mean(dat$coop[dat$likely.party=="R" & dat$state=="PA"],na.rm=T) + coef(m1)["likely.partyR"]*-1)*nrow(dat[dat$likely.party=="R" & dat$state=="PA" & dat$contact==1,])

#How many Independents in PA sample?
sum(dat$coop[dat$state=="PA" & dat$likely.party=="I"],na.rm=T)

#How many Independents if participation rate rose to level of Democrats?
(mean(dat$coop[dat$likely.party=="I" & dat$state=="PA"],na.rm=T) + coef(m1)["likely.partyR"]*-1)*nrow(dat[dat$likely.party=="I" & dat$state=="PA" & dat$contact==1,])



#######
#Figure2
#######

#Default Margins
par(mar=c(5.1,4.1,4.1,2.1))

state <- sort(unique(dat$state[!is.na(dat$likely.party)]))

rep.coop <- NA
ind.coop <- NA
rep.coop.se <- NA
ind.coop.se <- NA

for(i in 1:length(state)){
  m1<- felm(coop ~ likely.party + dem.perc16 + log(cases) +
              female +AgeUnder30  +Age3039 + Age4049 + Age5059 + Age6074 +
              white + black + hispanic + density, data=dat[dat$state==state[i],])
  rep.coop[i] <- coef(m1)["likely.partyR"]
  ind.coop[i]  <- coef(m1)["likely.partyI"]
  rep.coop.se[i]  <- sqrt(vcov(m1)["likely.partyR","likely.partyR"])
  ind.coop.se[i]  <-  sqrt(vcov(m1)["likely.partyI","likely.partyI"])
}

#View Results
cbind(state,rep.coop,ind.coop)

ypos1 <- seq(12,1,-1)
ypos2 <- seq(11.5,0.5,-1)
ypos3 <- seq(11.75,.75,-1)

pdf(file="./Figures/CLINTON-LAPINSKI-TRUSSLER-Figure2.pdf")
plot(rep.coop,ypos1, xlim=c(-0.20,0.20), ylim=c(0,13), axes=F, pch=16, col="firebrick",
     ylab="", xlab="Effect on P(Cooperate) vs. Democrat")
segments(rep.coop,ypos1, rep.coop +rep.coop.se*1.96,ypos1, col="firebrick")
segments(rep.coop,ypos1, rep.coop -rep.coop.se*1.96,ypos1, col="firebrick")


points(ind.coop, ypos2, col="purple", pch=17)
segments(ind.coop,ypos2, ind.coop +ind.coop.se*1.96,ypos2, col="purple")
segments(ind.coop,ypos2, ind.coop -ind.coop.se*1.96,ypos2, col="purple")


axis(side=1, at=seq(-0.20, 0.20,0.02))
axis(side=2, at=seq(11.75,.75,-1), labels=state, las=2)
abline(h=seq(1.25,11.25,1),lty=2, col="gray80")
abline(v=0, lty=3)
legend("topright", c("Republican","Independent/Other"), pch=c(16,17), col=c("firebrick","purple"))
dev.off()

#################
#Turnout Footnote
#################


mean(dat$turnout20,na.rm=T)
mean(dat$turnout20[dat$likely.party=="D"],na.rm=T)
mean(dat$turnout20[dat$likely.party=="R"],na.rm=T)
mean(dat$turnout20[dat$likely.party=="I"],na.rm=T)
table(is.na(dat$turnout20))




########
#Figure3
########



#States
state <- sort(unique(ep$state))
state <- rev(state)
margin.weighted1 <- rep(NA, length(state))
margin.weighted2<- rep(NA, length(state))
trump.actual <- rep(NA, length(state))
biden.actual <- rep(NA, length(state))

for(j in 1:length(state)){
  trump.actual[j] <- el$rep_percent[el$stateid==state[j]]
  biden.actual[j] <- el$dem_percent[el$stateid==state[j]]
}

trump.actual <- as.numeric(gsub("[\\%]","",trump.actual))
biden.actual <- as.numeric(gsub("[\\%]","",biden.actual))
margin.real <- biden.actual - trump.actual
survey.n <- rep(NA, length(state))

#Reweighting of polls, by state

for(j in 1:length(state)){
  #Generate state disposition data
  d <- dat[dat$state==state[j],]
  #Generate state exit poll data
  e <- ep[ep$state==state[j],]
  #Generate additional id, to play well with survey package
  e$id2 <- seq(1, nrow(e))
  #Save n of survey in state
  survey.n[j] <- nrow(e)

  #Use logit to get cooperation rates by party
  #(This returns the exact same result as calculating the means,
  #but using logit allows this to be scaled to more variables.)
  m <- glm(coop ~  likely.party, data=d, family = "binomial")

  #Make reduced dataset for prediction
  keep <- c("id2","likely.party")
  e.p<- e[keep]
  e.p <-e.p[complete.cases(e.p),]
  #Predict cooperation rate in the exit poll using logit result
  e.p$predicted.coop <- predict(m, type="response", newdata=e.p)

  #Invert predicted cooperation to get a cooperation weight
  e.p$coop.weight <- 1/(e.p$predicted.coop)

  #Merge this back into exitpoll data
  keep <- c("id2","coop.weight")
  e.p <- e.p[keep]
  e <- merge(e,e.p, by="id2")

  #Generate weighting objects
  original.weights <- svydesign(id=~id2, weights=~weight, data=e)
  e$reweight.coop <- e$weight*e$coop.weight
  coop.weights <- svydesign(id=~id2, weights=~reweight.coop, data=e)

  #Original Estimates
  trump <- as.numeric(prop.table(svytable(~vote.choice, design=original.weights)))[3]
  biden <- as.numeric(prop.table(svytable(~vote.choice, design=original.weights)))[1]
  margin.weighted1[j] <- biden - trump

  #Cooperation adjusted estimates
  trump <- as.numeric(prop.table(svytable(~vote.choice, design=coop.weights)))[3]
  biden <- as.numeric(prop.table(svytable(~vote.choice, design=coop.weights)))[1]
  margin.weighted2[j] <- biden - trump
}


margin.weighted1 <- round(margin.weighted1*100,2)
margin.weighted2 <- round(margin.weighted2*100,2)
margin.real <- round(margin.real,2)

ypos <- seq(1,length(state))

pdf(file="./Figures/CLINTON-LAPINSKI-TRUSSLER-Figure3.pdf")
plot(margin.real,ypos, xlim=c(-30,30), pch=16, xlab="Estimate of Biden - Trump Margin", ylab="",
     axes=F, col="black", type="p")
points(margin.weighted1, ypos, pch=17, col="darkorange")
points(margin.weighted2, ypos, pch=18, col="forestgreen")
segments(margin.weighted1, ypos,margin.weighted2, ypos, col="forestgreen", lwd=2)

abline(v=0, lty=2)
axis(side=1, at=seq(-30,30,10))
axis(side=2, at=ypos, labels=state, las=2)
legend("topright", c("Election Result","Original Survey Estimate", "Coop Reweighted Estimate"),
       pch=c(16,17,18), col=c("black","darkorange","forestgreen"), cex=.75)
dev.off()


improve <- (margin.weighted1 - margin.real) - (margin.weighted2 - margin.real)
weighted.mean(improve, w=1/survey.n)



#########
#Figure 4
#########


#A#
#Make figure that is proportion of Dems, Inds, and Reps that are correctly identified.

state <- sort(unique(ep$state))
state <- rev(state)

dem.correct <- rep(NA, length(state))
ind.correct <- rep(NA, length(state))
rep.correct <- rep(NA, length(state))

for(i in 1:length(state)){
  dem.correct[i] <-prop.table(table(ep$pid.lab[ep$state==state[i]], ep$likely.party[ep$state==state[i]]),2)[1,1]
  ind.correct[i] <-prop.table(table(ep$pid.lab[ep$state==state[i]], ep$likely.party[ep$state==state[i]]),2)[2,2]
  rep.correct[i] <-prop.table(table(ep$pid.lab[ep$state==state[i]], ep$likely.party[ep$state==state[i]]),2)[3,3]
}

cbind(state,dem.correct,ind.correct,rep.correct)

dem.correct <- dem.correct*100
ind.correct <- ind.correct*100
rep.correct <- rep.correct*100

ypos <- seq(1,length(state))
ypos2 <- ypos-.25
ypos3 <- ypos2-.25

pdf(file="./Figures/CLINTON-LAPINSKI-TRUSSLER-Figure4.pdf",
    height=6, width=8)
par(mfrow=c(1,2))
plot(dem.correct,ypos, xlim=c(0,100), ylim=c(0,7),pch=18, xlab="Percentage Likely Partisans Correctly Placed", ylab="",
     axes=F, col="dodgerblue", type="p")
points(ind.correct,ypos2, pch=17, col="purple", type="p")
points(rep.correct,ypos3, pch=16, col="firebrick", type="p")
abline(h=ypos-.8, lty=2)
abline(v=50, lty=3, lwd=2)
axis(side=1, at=seq(0,100,25))
axis(side=2, at=ypos2, labels=state, las=2)
legend("topleft", c("Democrat","Independent", "Republican"),
       pch=c(18,17,16), col=c("dodgerblue","purple","firebrick"), cex=.75)

#B
#Determine the degree to which likely partisanship is predictive of voe choice.
state <- sort(unique(ep$state))
state <- rev(state)

dem.biden <- rep(NA, length(state))
ind.biden <- rep(NA, length(state))
rep.biden <- rep(NA, length(state))

for(i in 1:length(state)){
  dem.biden[i] <-prop.table(table(ep$biden[ep$state==state[i]], ep$likely.party[ep$state==state[i]]),2)[2,1]
  ind.biden[i] <-prop.table(table(ep$biden[ep$state==state[i]], ep$likely.party[ep$state==state[i]]),2)[2,2]
  rep.biden[i] <-prop.table(table(ep$biden[ep$state==state[i]], ep$likely.party[ep$state==state[i]]),2)[2,3]
}

cbind(state,dem.biden,ind.biden,rep.biden)

dem.biden <- dem.biden*100
ind.biden <- ind.biden*100
rep.biden <- rep.biden*100

ypos <- seq(1,length(state))
ypos2 <- ypos-.25
ypos3 <- ypos2-.25

plot(dem.biden,ypos, xlim=c(0,100), ylim=c(0,7),pch=16, xlab="Percentage Likely Partisans Voting Biden", ylab="",
     axes=F, col="dodgerblue", type="p")
points(ind.biden,ypos2, pch=17, col="purple", type="p")
points(rep.biden,ypos3, pch=18, col="firebrick", type="p")
abline(h=ypos-.8, lty=2)
abline(v=50, lty=3, lwd=2)
axis(side=1, at=seq(0,100,25))
axis(side=2, at=ypos2, labels=state, las=2)

dev.off()













